Algorytm h-adaptacji
Celem algorytmu h adaptacji jest zwiększenie dokładności aproksymacji na siatkach obliczeniowych poprzez połamanie wybranych elementów na mniejsze. Algorytm h adaptacji możliwy jest do zaimplementowania dopiero wtedy kiedy skonstruujemy takie funkcje bazowe, które będzie można rozpinać na elementach skończonych różnego rozmiaru. W przypadku siatek niestrukturalnych, zbudowanych z elementów trójkątnych, możliwe jest zastosowanie algorytmu Rivary opisanego w innym rozdziale.
W przypadku siatek strukturalnych, zbudowanych z elementów prostokątnych, na których rozpięto funkcje bazowe zgodnie z zadanymi wektorami węzłów, h adaptacja jest możliwa wtedy jeśli pomiędzy wszystkimi łamanymi elementami wsadzimy \( C^0 \) separatory.
W rozdziale "Aproksymacja z wiszącymi węzłami" opisany został sposób łączenia ze sobą funkcji bazowych rozpiętych na sąsiadujących połamanych elementach różnego rozmiaru.
Algorytm h-adaptacji może zostać zaimplementowany w formie algorytmu a priori lub a posteriori, oraz stosując adaptacje izotropowe lub anizotropowe.
Algorytm adaptacyjny może używać jednej siatki obliczeniowej, tak zwanej siatki rzadkiej, lub dwóch siatek obliczeniowych rzadkiej i gęstej. W naszej klasyfikacji przez algorytm adaptacyjny a priori, rozumiemy algorytm który podejmuje decyzję przed skonstruowaniem i rozwiązaniem problemu na siatce gęstej, po rozwiązaniu problemu na aktualnej siatce rzadkiej. Z tego względu określenie "a priori" niejako "z góry" rozumiemy tutaj bez
dodatkowych obliczeń związanych na przykład z rozwiązywaniem problemu na dodatkowej siatce gęstej. Przez algorytm adaptacyjny a posteriori rozumiemy algorytm który podejmuje decyzje o adaptacje po wykonaniu pewnych dodatkowych obliczeń (nie tylko rozwiązaniu problemu na siatce rzadkiej) na przykład po wygenerowaniu siatki gęstej i rozwiązywaniu problemu na siatce gęstej.
- Algorytm h-adaptacji a priori używa jednej siatki obliczeniowej, a szacowanie błędu odbywa się na podstawie obliczeń lokalnych na poszczególnych elementach (na przykład możliwe jest szacowanie błedów na podstawie pierwszych lub/i drugich pochodnych kierunkowych rozwiązania lokalnie na elementach skończonych). Wiele metod wyboru adaptacji opisanych jest w pracy Marka Ainswortha i Johna Tinsleya Odena [1] (istnieje rozbieżność nazewnictwa metod pomiędzy naszym modułem a tym artykułem, wynikająca stąd iż u nas metodami a priory nazywane są metody używający jedynie siatki rzadkiej a metodami a posteriori metody używające dwóch siatek, u Autorów metoda używająca jednej siatki nazywa się a posteriori ponieważ jest ona uruchamiana po rozwiązaniu problemu, dla Autorów metoda a priori podejmowała by decyzje o adaptacji przed rozwiązaniem problemu obliczeniowego, patrząc na przykład na same współczynniki modelu i kształt siatki obliczeniowej).
- Algorytm h-adaptacji a posteriori używa dwóch siatek obliczeniowych, siatki rzadkiej i siatki gęstej, a szacowanie błędu na poszczególnych elementach odbywa się poprzez porównanie rozwiązań na siatkach rzadkiej i gęstej
- Algorytm h-adaptacji izotropowych, może działać w wersji a prior i i a posteriori, oraz dzieli on wybrane elementy na cztery (elementy prostokątne w dwóch wymiarach) lub osiem (elementy sześcienne w trzech wymiarach) nowe elementy (możliwe są też różne łamania elementów trójkątnych w dwóch wymiarach, oraz czworościennych, pryzmatycznych i piramid w trzech wymiarach; ostatnio zaczęto również używać metody elementów skończonych na wielobokach np. metoda PolyDPG)
- Algorytm h-adaptacji anizotropowych, może działać w wersji a priori i a posteriori, oraz dzieli on wybrane elementy na dwa nowe elementy w wybranym kierunku
Jak powiedział kiedyś prof. Zienkiewicz, jeden z twórców metody elementów skończonych, wystarczającą dokładność w wielu zagadnieniach inżynierskich dostaniemy stosując wielomiany kwadratowe i h-adaptację. Zakładamy więc że siatka początkowa to posklejane grupy elementów opisanych wektorami węzłów [0 0 0 1 1 1] i [0 0 0 1 1 1].
Algorytm h-adaptacji anioztropowej a priori w ogólnej postaci zapisać można w sposób następujący:
- niejednorodne h adaptacje wykonywane apropri (siatka początkowa, dokładność, maksymalna liczba kroków)
- aktualna siatka = siatka początkowa
- rozwiąż swój problem obliczeniowy na aktualnej siatce (na przykład projekcję bitmapy lub transport ciepła), obliczając aktualne rozwiązanie \( u \)
- pętla dla i = 1 do maksymalnej liczby kroków
- \( K_{maxerr }=0; K^x_{maxerr}=0; K^y_{maxerr}=0 \)
- dla każdego elementu \( K \) aktualnej siatki wykonaj następujące czynności
- oblicz błąd aproksymacji na elemencie \( K \) obliczając seminormę gradientu rozwiązania na elemencie \( K_{err}=\|\nabla(B - u) \|=\int |\frac{\partial (B-u)}{\partial x}|^2+ |\frac{\partial (B-u)}{\partial y}|^2dxdy \), oraz normy z pochodnych kierunkowych \( K^x_{err}=\|\nabla_x(B - u) \|=\int |\frac{\partial (B-u)}{\partial x}|^2dxdy \) oraz \( K^y_{err}=\|\nabla_y(B - u) \|=\int \|\frac{\partial (B-u)}{\partial y}|^2dxdy \)
- jeśli \( K_{err}>K_{maxerr} \) wtedy \( K_{maxerr}=K_{err} \)
- jeśli \( K^x_{err}>K^x_{maxerr} \) wtedy \( K^x_{maxerr}=K_{err} \)
- jeśli \( K^y_{err}>K^y_{maxerr} \) wtedy \( K^y_{maxerr}=K_{err} \)
- koniec pętli elementach
- jeśli \( K_{maxerr}> \) wymagana dokładność, to koniec
- pętla po elementach \( K \) siatki
- jeśli \( K_{err}>0.33*K_{maxerr} \) wówczas połam element w dwóch kierunkach, na cztery nowe elementy, i wygeneruj nowe lokalne funkcje bazowe, kontynuuj pętle po elementach
- jeśli \( K^x_{err}>0.33*K^x_{maxerr} \) wówczas połam element w kierunku osi \( x \) na dwa nowe elementy, i wygeneruj nowe funkcje bazowe, kontynuuj pętle po elementach
- jeśli \( K^y_{err}>0.33*K^y_{maxerr} \) wówczas połam element w kierunku osi \( y \) i wygeneruj nowe funkcje bazowe, kontynuuj pętle po elementach
- koniec pętli po elementach
- konieć pętli adaptacji
Poprzez błąd w przykładowym zadaniu projekcji bitampy, rozumiemy tutaj błąd w normie \( H^1 \) reprezentujący różnicę pomiędzy aproksymacją ciągłą np. odlewu terenu którego wysokość w różnych punktacj opisana jest za pomocą odcieni pikseli bitmapy.
Jeśli konieczne jest połamanie danego elementu (wykonanie h-adaptacji dla tego elementu), mamy do wyboru trzy scenariusze. Pierwsza możliwość to połamanie elementu na cztery mniejsze elementy. Druga możliwość to połamanie elementu na dwa nowe elementy w kierunku osi \( x \). Trzecia możliwość to połamanie danego elementu na dwa mniejsze elementy w kierunku osi \( y \).
W algorytmie adaptacyjnym zazwyczaj nie wykonuje się adaptacji na wszystkich elementach siatki, ale jedynie na tych elementach na których błąd jest większy niż 33 procent błedu maksymalnego (liczonego po wszystkich elementach). Jest to oczywiście arbitrarnie przyjęta wartość, dobrana na podstawie doświadczeń numerycznych przez prof. Leszka Demkowicza. W zaproponowanej wersji algorytmu próbujemy najpierw wykonać adaptacje w obydwu kierunkach (połamać element w obydwu kierunkach), a następnie (jeśli nie jest to wskazane) próbujemy połamać element w jednym z dwóch kierunków. Jest to oczywiście pzykładowa wersja algorytmu, możliwe są różne modyfikacje. Najbardziej zaawansowaną wersją algorytmu adaptacyjnego jest algorytm hp adaptacji opisany w innym module.
Generowanie nowych funkcji bazowych polega na generowaniu nowych lokalnych wektorów węzłów dla połamanych elementów. Na przykład, jeśli mamy element opisany wektorem węzłów [0 0 0 1 1 1] w kierunku osi \( x \) oraz wektorem węzłów [0 0 0 1 1 1] w kierunku osi \( y \), i połamiemy ten element na cztery nowe elementy, wówczas mamy wektor węzłów [0 0 0 0.5 0.5 1 1 1] w kierunku osi \( x \) oraz wektor węzłów [0 0 0 0.5 0.5 1 1 1] w kierunku osi \( y \).
Jeśli połamiemy ten element na dwa nowe elementy w kierunku osi \( x \), wówczas mamy wektor węzłów [0 0 0 0.5 0.5 1 1 1] w kierunku osi \( x \) oraz wektor węzłów [0 0 0 0 1 1 1] w kierunku osi \( y \).
Jeśli połamiemy ten element na dwa nowe elementy w kierunku osi \( y \), wówczas mamy wektor węzłów [0 0 0 1 1 1] w kierunku osi \( x \) oraz wektor węzłów [0 0 0 0.5 0.5 1 1 1] w kierunku osi \( y \).
Scalając ze sobą wiele patchów (grup elementów) na których rozpięte są różne wielomiany B-spline, uzyskamy siatkę obliczeniową na której możliwe będą lokalne modyfikacje stopnia funkcji B-spline, na przykład stosując opisany algorytm p-adaptacji. Sposób scalania funkcji bazowych uzyskanych na dużych i małych połamanych elementach został opisany w rozdziale "Aproksymacja z wiszącymi węzłami".
Algorytm automatycznej h-adaptacji aposteriori:
- Algorytm automatycznej h adaptacji (siatka początkowa, wymagana dokładność, maksymalna liczba iteracji)
- siatka rzadka = siatka początkowa
- pętla i = 1 do maksymalnej liczby iteracji
- Rozwiąż problem obliczeniowy na aktualnej siatce rzadkiej (na przykład problem projekcji bitmapy lub problem transportu ciepła) uzyskując rozwiązanie \( u_h \)
- siatka gęsta = siatka rzadka
- Połam każdy element siatki rzadkiej na cztery elementy
- Rozwiąż problem obliczeniowy na aktualnej siatce gęstej uzyskując rozwiązanie \( u_{h/2} \)
- Bład maksymalny \( K_{max}=0 \)
- Pętla po elementach \( K \) siatki rzadkiej
- Dla każdego elementu siatki rzadkiej oszacuj błąd względny (normę różnicy pomiędzy rozwiązaniem na siatce rzadkiej a gęstej) \( K_{rel}=\| u_h - u_{h/2} \|_K = \int_K \left( u_h-u_{h/2} \right)^2 + \left( \frac{\partial u_h-u_{h/2}}{\partial x} \right)^2 + \left( \frac{\partial u_h-u_{h/2}}{\partial y} \right)^2 \)
- Jeśli \( K_{rel}>K_{max} \) wówczas \( K_{max}=K_{rel} \)
- Koniec pętli po elementach
- Jeśli \( K_{max}> \) wymaganej dokładności, wówczas koniec.
- Pętla po elementach \( K \) siatki rzadkiej
- Jeśli \( K_{rel}>0.33 K_{max} \) wówczas wybierz optymalny sposób adaptacji elementu z siatki rzadkiej i zaaplikuj go do elementu \( K \) siatki rzadkiej
- Koniec pętli po elementach
- Koniec iteracji
Mamy więc trzy możliwości adaptacji pojedynczego elementu:
- Złamanie elementu na cztery nowe elementy
- Złamanie elementu na dwa nowe elementy w kierunku poziomym
- Złamanie elementu na dwa nowe elementy w kierunku pionowym
W jaki sposób podejmowana jest decyzja o wyborze rodzaju adaptacji pojedynczego elementu?
Decyzja podejmowana jest zgodnie z następującym algorytmem.
- Wybór optymalnej strategii adaptacji elementu \( K \) (rozwiązanie na siatce rzadkiej zawężone do elementu \( u_h \), rozwiązanie na siatce gęstej zawężone do elementu \( u_{h/2} \))
- Pętla po możliwych rodzajach adaptacji elementu od 1 do 3
- Największa prędkość spadku błedu = 0
- Optymalna adaptacja = 0
- Wykonaj kopię \( J \) elementu \( K \) i wykonaj na niej rozważaną adaptację
- Oblicz projekcję \( w \) rozwiązania na siatce gęstej \( u_{h/2} \) na adaptowany element \( J \)
- Oblicz o ile spadnie błąd względny jeśli wykonamy adaptację reprezentowaną przez kod, spadek błędu (kod)= \( \|u_{h/2}-u_h\|_K-\|u_{h/2}-w\|_K \)
- Oblicz ile niewiadomych (ile funkcji bazowych) musimy dodać żeby zrealizować adaptację reprezentowaną przez kod, koszt adaptacji (kod)
- Oblicz i zapamiętaj prędkość spadku błędu (kod) = spadek błędu (kod) / koszt adaptacji (kod)
- Jeśli prędkość spadku błędu (kod) jest większa niż największa prędkość spadku błedu, to zapamiętaj największą wartość spadku błędu = pedkość spadku błędu (kod), optymalna adaptacja=kod
- Koniec pętli po możliwych rodzajach adaptacji
- Wykonaj na elemencie \( K \) znalezioną optymalną adaptację
Innymi słowy, wybieramy taki rodzaj adaptacji dla elementu, który pozwala uzyskać największy spadek błędu najmniejszym kosztem. Wielkość ta reprezentowana jest przez prędkość spadku błędu, która rośnie wraz ze spadkiem błędu, ale maleje wraz z poniesionym kosztem (z liczbą funkcji dodaną do elemetu, ponieważ koszt obliczenia na danym elemencie zależy od liczby niewiadomych, które są współczynnikami funkcji bazowych).